function [xs,apstar,cstar, aptilde, ctilde] = bc_opt_pol_iid_egm(bta,r,b,lybar,sigma_ly,gam, amin, amax, na, util_type,p_unemp, w_unemp)

%  iteration parameters
tol = 1e-4;
maxiter = 1e+4;

%  fix exogenous state y grid
ny = 101;
m = 6;
lyrange = m*sigma_ly;
lys = linspace(lybar-lyrange, lybar+lyrange, ny)';

%  fix end-of-period state a' grid
as = linspace(amin, amax, na)';

%  probabilities over log(y') (iid draws from N(lybar,sigma_y^2))  
%  TEST WHETHER QUADRATURE METHOD IS FASTER

w = lys(2) - lys(1);
prob = normcdf(lys + w/2, lybar, sigma_ly) - normcdf(lys - w/2, lybar, sigma_ly);
prob(1) = normcdf(lys(1) + w/2, lybar, sigma_ly);
prob(end) = 1 - normcdf(lys(end) - w/2, lybar, sigma_ly);
ys = exp(lys);

if p_unemp > 0
    ys = [ w_unemp*sum(ys.*prob); ys]; 
    prob = [p_unemp; (1-p_unemp)*prob]; 
    ny = length(ys); 
end

%  iterations

i = 1;
conv = 1;
convtilde = 1;
%ys = [w_unemp
[ass, yss] = ndgrid(as,ys);
ap = ass;   % initial guess of a'(a,y) at gridpoints (as,ys)
aptilde = ap;

while i < maxiter && conv > tol
    
    %  consumption from Euler eqn
    if strcmp(util_type,'crra')
        val = (ass.*(1+r) + yss - ap).^(-gam); % crra utility
        c = (bta*(1+r)*val*prob).^(-1/gam); % crra utility
    elseif strcmp(util_type, 'quad')
        val = 1-gam*(ass.*(1+r) + yss - ap); % quad utility
        c = (1-bta.*(1+r).*val*prob)./gam; % quad utility
    else
        error('unknown utility type');
    end
    
    % this is the end. grid for a implied by the end. grid of c
    
    a_end = (ass - yss + repmat(c, 1,ny))./(1+r);
    
    %  update policy function a'(a,y) cycling through values of y
    ap_new = NaN(size(ass));
    aptilde_new = NaN(size(ass));
    for j=1:ny
        
        % find policy function on (y',a') by interpolation
        ap_new_temp  = interp1(a_end(:,j), ass(:,j), ass(:,j), 'linear','extrap');
        ap_new(:,j) = ap_new_temp;
        
        % unconstrained optimal action today, given tomorrow's (constrained) policy
        aptilde_new(:,j) = ap_new(:,j);
        
        % if policy function violates borrowing constraint, make it
        % bind
        ap_new_temp(ap_new_temp<b) = b;
        ap_new(:,j) = ap_new_temp;
        
    end
    
    %  convergence
    
    conv = max(max(abs((ap_new - ap)./(0.5.*(ap_new+ap)))));
    convtilde = max(max(abs(aptilde_new - aptilde)));
    i = i + 1;
    ap = ap_new;
    aptilde = aptilde_new;
    
    %  print iteration number
    %if mod(i,100)==0
    %    fprintf('\t egm recursion = %d\t conv = %d\n', i, conv)
    %end
end
%fprintf('\t egm recursion = %d\t conv = %d\n', i, conv)

%  policy functions

apstar = ap(:);
aptilde = aptilde(:);
xs = ass(:).*(1+r) + yss(:);
cstar = xs - apstar;
ctilde = xs - aptilde;
end
